Skip to contents

happign for foresters

Source: vignettes/web_only/happign_for_foresters.Rmd
happign_for_foresters.Rmd

First choose a zone of interest

For the example we will work with the Gâvre forest. The starting point come from google map : you just have to select a random point inside the forest by using right click on map.

point <- st_sfc(st_point(c(-1.7977020162748798, 47.52634370369097)), crs = 4326)

tm_shape(point)+
   tm_dots(size = 0.1)+
   tm_basemap("OpenStreetMap")

Extract commune borders

Now that we have point, getting shape of commune is just a matter of finding the good layer name from IGN. The one we are after is called “ADMINEXPRESS-COG-CARTO.LATEST:commune” from the “Administrative” category. Then we use the get_wfs() function to download the shape.

As you can see below, there are two communes close to our point (tip : click on the interactive map below to get more info). For following examples we will only use “Le Gâvre” shape.

apikey <- get_apikeys()[1]  # "administratif"
layer_name <- "ADMINEXPRESS-COG-CARTO.LATEST:commune"

# The layer_name I use come from `all_layers = get_layers_metadata(apikey, "wfs") that return all of them layer_name

borders <- get_wfs(shape = point,
                   apikey = apikey,
                   layer_name = layer_name)
#> 1/1 downloaded

tm_shape(borders)+
   tm_polygons(lwd = 2,
               col = "nom",
               id = "nom",
               popup.vars = "population")+
   tm_basemap("OpenStreetMap")

Download cadastral parcel

Cadastral parcels are essential for any forest management. Here how to download it with happign.

# Select only commune to reduce downloading time for hight resolution raster
borders <- borders %>% 
   slice(1)

apikey <- get_apikeys()[14] # "parcellaire"
layer_name <- get_layers_metadata(apikey, "wfs")
name_parcellaire_layer <- as.character(layer_name[15,2]) # "CADASTRALPARCELS.PARCELLAIRE_EXPRESS:parcelle"

parcellaire <- get_wfs(shape = borders,
                       apikey = apikey,
                       layer_name = name_parcellaire_layer) %>% 
   st_transform(st_crs(borders)) %>% 
   st_intersection(borders)
#> 1/10 downloaded
#> 2/10 downloaded
#> 3/10 downloaded
#> 4/10 downloaded
#> 5/10 downloaded
#> 6/10 downloaded
#> 7/10 downloaded
#> 8/10 downloaded
#> 9/10 downloaded
#> 10/10 downloaded
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

tm_shape(borders)+
   tm_borders(col = "red", lwd = 2)+
tm_shape(parcellaire)+
   tm_polygons(alpha = 0)+
   tm_basemap("OpenStreetMap")

Rq :

  • IGN WFS can return a maximum of 1000 features. The get_wfs() function overrides this limit by performing several consecutive requests, that why there 10 request aboves ;
  • The shape used get_wfs() is converted to bbox before requesting. In other words, several parcels will be outside our zone, which is why an intersection is made.

Downloading BD Forêt

The first interesting layer for forester is the “BD Forêt” which is all vegetation type assigned to each area greater than or equal to 0.5 ha (5,000 m²). There is two layer for forest : the old one BD Forêt V1 and the new one BD Forêt V2.

apikey <- get_apikeys()[9]  #"environement"
layer_name <- get_layers_metadata(apikey, "wfs")
name_BDV1 <- as.character(layer_name[1,2]) #LANDCOVER.FORESTINVENTORY.V1:resu_bdv1_shape
name_BDV2 <- as.character(layer_name[2,2]) #LANDCOVER.FORESTINVENTORY.V2:bdforetv2

BDV1 <- get_wfs(borders, apikey, name_BDV1) %>% 
   st_transform(st_crs(borders)) %>% 
   st_intersection(borders) # Layer is download for a bbox around the shape. An intersection is done to get forest only inside our zone of interest
#> 1/1 downloaded

BDV2 <- get_wfs(borders, apikey, name_BDV2) %>%
   st_transform(st_crs(borders)) %>% 
   st_intersection(borders)
#> 1/1 downloaded

tm_shape(BDV1) +
   tm_polygons(col = "libelle",
               popup.vars = names(BDV1)[1:20],
               legend.show = FALSE) +
tm_shape(BDV2) +
   tm_polygons(col = "tfv",
               alpha = 0.5,
               popup.vars = names(BDV2)[1:16],
               legend.show = FALSE) +
tm_shape(borders) +
   tm_borders(lwd = 2) +
   tm_basemap("OpenStreetMap")

With the map below, you can choose which BD Forêt you want to display. The BDV2 is obviously more precise. Also, many differences can be calculated, one is the area : BDV2 as 285737.31940975 more hectare described.

More calculations can be done as you can see below :

forest_type_BDV2 = BDV2 %>%
  mutate(area = as.numeric(st_area(geometry))) %>%
  st_drop_geometry() %>%
  group_by(essence) %>%
  summarise(sum_area = sum(area)/10000) %>%
  arrange(sum_area) %>%
  mutate(essence = as.factor(essence))

ggplot()+
  geom_col(data = forest_type_BDV2,
           aes(x = rev(reorder(essence, sum_area)),
               y = sum_area,
               fill = as.factor(essence)))+
  theme_bw()+
  labs(title = "Surface couverte par essences [ha]",
       y = "Surface [ha]",
       fill = "Essence :")+
  theme(axis.text.x = element_blank())

Detect protected area

One information you really want when you work at forest management is if your zone of interest is inside protected area. The example code below is design to automatically test every layer starting with “PROTECTED” so you can be sure that you have all of them.

So, for Le Gâcre forest we find : - 19 forest rescue point (PRSF in red) - 1 SIC (“Sites d’Importance Communautaire” in purple) - 1 ZPS (“Zone de protection spéciale” in green) - 6 Znieff 1 (blue on map) - 3 Znieff 2 (grey on map)

Again, you can click on map, point and shape for more informations.

apikey <- "environnement"
protected_area_names <- get_layers_metadata(apikey, "wfs") %>% 
   filter(grepl("^PROTECTED", name)) %>% 
   pull(name)

all_protected_area <- map(.x = protected_area_names,
                         .f = ~ try(get_wfs(borders, apikey, .x), silent = TRUE)) %>% 
   set_names(protected_area_names) %>% 
   discard(~ nrow(.) == 0)

# Plot the result
tm_shape(all_protected_area[[1]])+
   tm_dots(group = "Point rencontre des secours en forêts", col = "red")+
tm_shape(all_protected_area[[2]])+
   tm_polygons(group = "SIC", alpha = 0.8, col = "purple")+
tm_shape(all_protected_area[[4]])+
   tm_polygons(group = "Znieff 2", alpha = 0.8)+
tm_shape(all_protected_area[[3]], alpha = 0.8)+
   tm_polygons(group = "Znieff 1", alpha = 0.8, col = "blue")+
tm_shape(all_protected_area[[5]])+
   tm_polygons(group = "ZPS", alpha = 0.8, col = "green")+
tm_shape(borders) +
   tm_borders(lwd = 2) +
   tm_basemap("OpenStreetMap")

MNS, MNT and MNH…

It’s always good to know a more about terrain topologie. IGN offers MNT and MNS for download. As a reminder, the MNT corresponds to the surface of the ground and the MNS to the real surface (in our case, the trees). It is thus easy to find the height of the trees by subtracting the DTM from the MNS.

layers_name <- get_layers_metadata("altimetrie", "wms")
MNT_name <- layers_name[2,2] # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES"
MNS_name <- layers_name[3,2] # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES.MNS"

MNT <- get_wms_raster(borders, "altimetrie", resolution = 5, MNT_name, filename = "MNT")
#> 1/2 downloaded
#> 2/2 downloaded
file.remove("MNT_5m.tif") # I don't want to keep MNT on disk because this is vignette
#> [1] TRUE
MNS <- get_wms_raster(borders, "altimetrie", resolution = 5, MNS_name, filename = "MNS")
#> 1/2 downloaded
#> 2/2 downloaded
file.remove("MNS_5m.tif") # I don't want to keep MNS on disk because this is vignette
#> [1] TRUE

MNH <- MNS - MNT
MNH[MNH < 0] <- NA  # Remove negative value 
MNH[MNH > 50] <- NA # Remove height more than 50m

tm_shape(MNH)+
   tm_raster()+
tm_shape(borders)+
   tm_borders(lwd = 2)+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 1193 by 839 cells. See tm_shape manual (argument raster.downsample)

NDVI

The code below present the calculation of the NDVI. All informations and palette come from this website The value range of an NDVI is -1 to 1. It is (Near Infrared - Red) / (Near Infrared + Red) :

  • Water has a low reflectance in red, but almost no NIR (near infrared) reflectance. So the difference will be small and negative, and the sum will be small, and NDVI large and negative.
  • Plants have a low reflectance in red, and a strong NIR reflectance. So the difference will be large and positive, and the sum will be just about the same as the difference, so NDVI will be large and positive.

Categories are somewhat arbitrary, and you can find various rules of thumb, such as:

  • Negative values of NDVI (values approaching -1) correspond to water. Values close to zero (-0.1 to 0.1) generally correspond to barren areas of rock, sand, or snow. Low, positive values represent shrub and grassland (approximately 0.2 to 0.4), while high values indicate temperate and tropical rainforests (values approaching 1).
  • Very low values of NDVI (0.1 and below) correspond to water, barren areas of rock, sand, or snow. Moderate values represent shrub and grassland (0.2 to 0.3), while high values indicate temperate and tropical rainforests (0.6 to 0.8).
# To show the possibility of 20cm resolution for IRC, let's take only the biggest parcels
biggest_parcels <- parcellaire %>% 
   mutate(area = st_area(.)) %>% 
   slice_max(area)

IRC <- get_wms_raster(shape = biggest_parcels,
                      apikey = "ortho",
                      resolution = 0.2,
                      layer_name = "ORTHOIMAGERY.ORTHOPHOTOS.IRC",
                      filename = "IRC")
#> 1/4 downloaded
#> 2/4 downloaded
#> 3/4 downloaded
#> 4/4 downloaded
file.remove("IRC_0_2m.tif") # I don't want to keep IRC on disk because this is vignette
#> [1] TRUE

infrared <- IRC %>% slice(band, 1)
red <- IRC %>% slice(band, 2)

NDVI <- (infrared-red)/(infrared+red)
   
breaks_NDVI <- c(-1,-0.2,-0.1,0,0.025 ,0.05,0.075,0.1,0.125,0.15,0.175,0.2 ,0.25 ,0.3 ,0.35,0.4,0.45,0.5,0.55,0.6,1)
   
palette_NDVI <- c("#BFBFBF","#DBDBDB","#FFFFE0","#FFFACC","#EDE8B5","#DED99C","#CCC782","#BDB86B","#B0C261","#A3CC59","#91BF52","#80B347","#70A340","#619636","#4F8A2E","#407D24","#306E1C","#216112","#0F540A","#004500")

tm_shape(borders)+
   tm_borders(lwd = 2, col = "red")+
tm_shape(NDVI)+
   tm_raster(stretch.palette = F,
             breaks = breaks_NDVI,
             palette = palette_NDVI,
             colorNA = "red")+
tm_shape(biggest_parcels)+
   tm_borders(lwd = 2, col = "blue")+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 1255 by 797 cells. See tm_shape manual (argument raster.downsample)

The gloss index

The gloss index represents the average of the image glosses. This index is therefore sensitive to the brightness of the soil, related to its moisture and the presence of salts on the surface. It characterizes especially the albedo (solar radiation that is reflected back to the atmosphere). The gloss index allows us to estimate whether the observed surface feature is light or dark.


gloss_index <- sqrt(red^2 + infrared^2)

tm_shape(borders)+
   tm_borders(lwd = 2, col = "red")+
tm_shape(gloss_index)+
   tm_raster()+
tm_shape(biggest_parcels)+
   tm_borders(lwd = 2, col = "blue")+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 1255 by 797 cells. See tm_shape manual (argument raster.downsample)

Last but not least… BD Topo

BD topo from IGN covers in a coherent way all the geographical and administrative entities of the national territory. So you can find in it :

  • Administrative (boundaries and administrative units);
  • Addresses (mailing addresses) ;
  • Building (constructions) ;
  • Hydrography (water-related features) ;
  • Named places (place or locality with a toponym describing a natural space or inhabited place);
  • Land use (vegetation, foreshore, hedge);
  • Services and activities (utilities, energy storage and transportation, industrial sites);
  • Transportation (road, rail or air infrastructure, routes);
  • Regulated areas (most of the areas are subject to specific regulations).

For the example below I choose to download all water-related data :

cour_eau <- get_wfs(borders, "topographie", "BDTOPO_V3:cours_d_eau") %>% 
   st_intersection(st_buffer(borders,100))
#> 1/1 downloaded
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
detail_hydro <- get_wfs(borders, "topographie", "BDTOPO_V3:detail_hydrographique") %>% 
   st_intersection(st_buffer(borders,100))
#> 1/1 downloaded
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# water detected by satellite
surf_hydro <- get_wfs(borders, "topographie", "BDTOPO_V3:surface_hydrographique") %>% 
   st_intersection(st_buffer(borders,100))
#> 1/1 downloaded
#> st_as_s2(): dropping Z and/or M coordinate
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

tm_shape(cour_eau)+
   tm_lines(col = "blue")+
tm_shape(detail_hydro)+
   tm_dots(col = "red")+
tm_shape(surf_hydro)+
   tm_polygons("steelblue")+
tm_shape(borders)+
   tm_borders(lwd = 2)+
   tm_basemap("OpenStreetMap")

What about history ?

The “Etat-major” map is a general map of France made, in its first version, in the 19th century. Here how to get it :

etat_major <- get_wms_raster(shape = borders,
                            apikey = "cartes",
                            resolution = 5,
                            layer_name = "GEOGRAPHICALGRIDSYSTEMS.ETATMAJOR40",
                            filename = "etat_major")
#> 1/2 downloaded
#> 2/2 downloaded
file.remove("etat_major_5m.tif") # I don't want to keep etat_major on disk because this is vignette
#> [1] TRUE

tm_shape(etat_major)+
   tm_rgb()+
tm_shape(borders)+
   tm_borders(lwd = 3, col = "red")+
   tm_basemap("OpenStreetMap")
#> stars object downsampled to 1193 by 839 cells. See tm_shape manual (argument raster.downsample)